The dynamics of dark solitons in a trapped superfiuid Fermi gas 
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We study soliton oscillations in a trapped superfiuid Fermi gas across the Bose-Einstein con- 
densate to Bardeen-Cooper-Schrieflfer (BEC-BCS) crossover. We derive an exact equation for the 
oscillation period in terms of observable quantities, which we confirm by solving the time-dependent 
Bogoliubov-de Gennes equations. Hence we reveal the appearance and dynamics of solitons across 
the crossover, and show that the period dramatically increases as the soliton becomes shallower on 
the BCS side of the resonance. Finally, we propose an experimental protocol to test our predictions. 

PACS numbers: 03.75.Ss, 03.75.Lm 



Solitons have been the focus of much recent research in 
the field of cold atoms, due to their ubiquitous production 
in the dynamics of BECs [1| . Their different forms cre- 
ate a broad family, from the common "grey" and "black" 
solitons in repulsive BECs, to the "bright" solitons in 
attractive BECs and "gap" solitons in optical lattices, 
and their more exotic cousins such as the "bright-dark" 
solitons, which were recently observed in two-component 
BECs We expect solitons to play an equally impor- 

tant role in the dynamics of degenerate Fermi gases. Even 
more fundamentally, topological excitations offer insights 
into the nature of coherence and superfluidity across the 
BEC-BCS crossover, as illustrated by the recent observa- 
tion of vortex lattices Despite this interest, the na- 
ture of soliton dynamics in Fermi gases remains elusive, 
and only stationary "black" solitons have been simulated 
across the BEC-BCS crossover Q. 

In this Letter, we investigate the oscillation of soli- 
tons in a harmonic trap. From fundamental statements 
about the nature of the soliton and the media it moves in, 
we derive universal relations, valid for both Bosonic and 
Fermionic superfluids, relating the soliton energy and os- 
cillation period T s to observable quantities such as phase 
jump, speed and density. The special case of unitarity 
is particularly interesting because the soliton width is of 
the order of interatomic distances and its observation will 
give access to the short-range physics. We then perform 
numerical simulations of the time-dependent Bogoliubov- 
de Gennes (TDBdG) equations 0] ■ By extracting the ap- 
propriate observable quantities from our simulations we 
find good agreement between the numerical and analytic 
models, which show that T s increases dramatically as the 
soliton becomes shallower when moving from a BEC to 
a BCS regime. Finally, we propose and simulate an ex- 
perimental protocol to demonstrate the variation in T s 
across the BEC-BCS crossover. 

General theory. Let us consider a soliton in a super- 
fluid gas, which may be either Bosonic or Fermionic, con- 
fined in an elongated trap with axial potential U (x) — 
mWji 2 /2, in which m is the atomic mass and uj x is the an- 
gular frequency. We assume that the width of the cloud 
in the axial direction is large compared to the size of 



the soliton. Then, in a local density approximation, the 
soliton can be treated as a point-like particle at position 
X\ its dynamics can be formulated in terms of the soliton 
energy in a uniform fluid E s (p, V 2 ) , where fi is the chem- 
ical potential of the fluid and V — dX/dt is the soliton 
velocity. Furthermore, we may say that \x (X) = p (0) — 
U (X). We define the inertial mass mj = 2dE s /dV 2 and 
the number N s — —dE s /dp — J_ ao [nid(x) — nido] dx, in 
which nidix) is the one-dimensional density and nido is 
the bulk density far from the soliton. The quantity N s 
is the deficit of particles associated with the depression 
in density at the soliton position. Note that typically 
mi < and N s < 0. Then energy conservation in the 
absence of dissipation gives @, H| 

dE s _ dE s dp (X) dX dE s dV 2 dV _ 

~dt ~ dp (X) dX dT + W 2 ~dV~dt~ ^' 

and hence m I (dV/dt) = -N s (dU/dX) = -N s muj 2 x X. 
For small amplitude oscillations, the soliton period is 

T s = y/rm/NsmTx, (2) 

where T x = 2tt/u> x , and mj and N s are taken for V — > 0. 

The key to our analytic treatment of solitons is to 
recognise the distinction between the "physical" momen- 
tum of the soliton P s — m J_ oc jdx, associated with the 
current 3 carried by the wave function, and the canonical 
momentum of the soliton P c . By performing a Galilean 
transformation into the frame of the soliton moving at 
velocity V, we find the current in the soliton frame 
3=3— n\dV ■ Since 3 = — n\doV , we derive 

/CO 
[n ld (x) - moo] dx = mVN s . (3) 
-00 

The canonical momentum of the soliton is instead defined 
by V = 8E s /dP c . It follows that 

dP c 1 dE s BE S 
The momenta P s and P c are different because, despite 



being a localized object from the point of view of the 
density profile and the velocity field, the soliton creates 
a jump J v in the phase ip of the order parameter. This 
phase jump is exploited when creating solitons in ex- 
periment with the "phase imprinting" technique 0, [t^- 
1 1 XI ] - In any real experiment, where the soliton is created 
by a localized perturbation, this phase difference will be 
compensated by a "counterflow" , which carries an addi- 
tional momentum AP. The difference between P s and 
P c was first found in a BEC described by the Gross- 
Pitaevskii (GP) equation and its meaning was dis- 
cussed in Ref . [l3[ • Far from the soliton we may say that 
the velocity field is v = fiVtp/mB, in which ms = m 
for Bosons and raj = 2m for Fermions. Hence we find 
AP = nidom f vdx — KnidoJ v m/mB, and [3] 



P c = P s + AP = P s + hn ldQ J v m/m B . 



(5) 



Taking into account that V — for J v — tt, and using 
Eqs. ([3j> and (J4j) , we obtain the important relation 



mN,V - 2 



BE, 
dV 2 



dV = —Hmdo(J<p — n)m/mB ■ (6) 



We then differentiate both sides of Eq. (|6]) with respect 
to V to derive 



d (N B V) hni do m dJ v 

m mi = — 

dV m B dV 



(7) 



Substituting m; using Eq. (|2J), for V — > we obtain the 
final result for T s : 



tJ ~ m B N s dV 



(8) 



Note that Eq. (jH)) contains only quantities which can 
be directly measured in numerical and real experiments. 
We stress that the basic assumptions used to derive this 
result are the absence of dissipation and the existence 
of an order parameter obeying Galilean invariance for 
particles of mass tub ■ 

Bogoliubov - de Gennes simulations. To test Eq. ([5]l. 
we model the dynamics of a three-dimensional (3D) su- 
perfluid Fermi gas by solving the TDBdG equations @ 



H 


A(r,t)~ 




u v (r,t) 




u v (r,t) 


a*M) 


-H 




_ v v (r,t) _ 


= lh dt 





(9) 

where H = —ti 2 V 2 /2m + U—fj,(0) and the order param- 
eter A is calculated as A = —gJ2 v u v v ri> m wn i cn 9 i s 
given by l/k f a = 8irE f /gk 3 f + ^AE c /it 2 Ef Here 
a is the 3D s-wave scattering length characterizing the 
interaction between atoms of different spins, while Et = 

1 /3 

h 2 k 2 /2m and kf = (37r 2 n) are the Fermi energy and 
momentum of an ideal Fermi gas of density n. The cut-off 
energy E c is introduced in order to remove the ultraviolet 
divergences in the TDBdG equations with contact poten- 
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FIG. 1: (a): Soliton oscillating in the density profile of a 
40 K superfluid for l/kfa = —0.5 with uj x = 27r50 rad s" 1 , 



3.3 ^im and a peak density n p 
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Enlargement of region contained within the dashed white box 
in (a), (c) & (d): Corresponding enlargements for l/kfa = 
0.0 and 1.0. Left and right insets in (b), (c) and (d) show 
the real part |5R (A)| and imaginary part |SJ (A)| of the order 
parameter in the regions contained within the dotted white 
boxes. Horizontal bars show scale. 



tials. The density of the gas is n(r,t) — 2 XL |^(r,t)| . 
Since the potential U has no y or z dependence, we write 
the BdG eigenfunctions as u v (r) = u ri {x)e l ' Kkyy+kzZ ^ and 
iV ( (r) = v ri {x)e l ( kyy+kzZ \ in which k y and k z are quan- 
tized according to k y = 2Tra y /L± and k z = 2ira z /L± 1 
where a y and a z are integers and L± is the width of the 
box in the y- and z-directions. 

As initial states at t = 0, we find stationary solutions 
of Eq. © • This has been done to investigate stationary 
black solitons across the BEC-BCS crossover [f| . We use 
the same technique to generate momentarily stationary 
solitons away from the trap center. When such an initial 
state is evolved in time, the black soliton is accelerated 
by the trap potential, and becomes grey. 

Figure Q] presents three typical simulations of the 
TDBdG equations. Panel (a) shows a soliton oscillat- 
ing in the density profile of a 40 K superfluid for l/kfa = 
—0.5, with lo x = 27r50 rad s _1 , L± = 3.3 /xm and a peak 
density n p = 1.8 x 10 18 m~ 3 = n ld {Q)/L\. Panel (b) 
is an enlargement of the central region of panel (a) con- 
tained within the dashed white box. Panels (c) and (d) 
are equivalent enlargements for l/kfa = and 1 respec- 
tively. We take E c = 30E f (E c = 50E f ) for l/kfa < 
{l/kfa > 0). For l/k f a = -0.5 [Figs. Ufa) and (b)], 
the soliton creates a shallow depression in the density 
of the cloud, on either side of which are smaller oscil- 
lations, known as Friedel oscillations. We also plot the 
real part |5R(A)| and imaginary part |9f(A)| of the or- 
der parameter in the region contained within the dotted 
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FIG. 2: T s (a) and N s (b) plotted against 1/fc/o. We find that 
T a = 1.7Ta; at unitarity. White (light gray, dark gray, black) 
crosses denote data for E c = 30% (50%, 75%, 100%). N s 
is calculated for transverse width Lj_ = 3.3 fim. 



white box in Fig. QJ b) as left and right insets respectively. 
Initially |3(A)| is zero, indicating that J v — n. As the 
soliton accelerates, the density depression becomes shal- 
lower, |3(A)| increases from zero and J v reduces. The 
insets also show that both |Sft(A)| and |9f(A)| contain 
Friedel oscillations in the vicinity of the soliton [17j . as 
in the density profile. This is in contrast to solitonic so- 
lutions of the GP equation, which always have a constant 
imaginary component of the order parameter [l8j . 

As l/kfa increases, the profile of the solitons in n and 
A tends to that of GP solitons. At unitarity [Fig. [He)], 
the Friedel oscillations are barely visible in n, but |3 (A)| 
(right inset) still contains a small dip at the position of 
the soliton. When l/k f a reaches 1.0 [Fig.^d)], |3(A)| 
(right inset) is almost constant across the cloud. We also 
observe that the density depression becomes deeper. For 
l/kfa = 1.0, the density minimum in the soliton is close 
to zero when it is stationary at the apex of an oscillation. 

Figure [T] also illustrates that the period T s decreases 
as we move from the BCS to the BEC regime. This 
effect is quantified in Fig. [2][a), which plots T s against 
l/kfa. The graph shows that T s drops rapidly as we 
move from the BCS to unitary regimes, before tending to 
the GP prediction of \f2T x in the BEC limit of large 
l/kfa. Our prediction of T s = (1.7±0.05)Tj; at unitarity 
is consistent with the value of V3T X given in Ref. fl7j . 
It is computationally difficult to reach convergence for 
large l/kfa, because a large number of states must be 
included in order to describe the formation of Bosonic 
molecules. To illustrate the gradual convergence towards 
the GP prediction, we plot three points for l/kfa = 1 
with E c = 50%, 75% and 100%, with a light gray, 
dark gray and black cross respectively. We also note that 
T s for l/kfa = —0.5 is lower than expected by looking 
at the general trend. This is a real effect that occurs 
because the pair size is becoming comparable with the 
width of the cloud. It can be avoided by reducing uj x . 

To compare the numerical results with the prediction 
of Eq. ©, we must also calculate N s and dJ v /dV. The 
quantity N s may be determined from stationary solu- 
tions of Eq. ©• We plot results in Fig. Hfb) as a func- 




FIG. 3: [{T s /T x f - 1] plotted against 1/ |7V S |. White (light 
gray, dark gray, black) symbols denote data for E c = 30% 
(50%, 75%, 100%). Solid line shows the prediction of 
Eq. (|11|) . Inset shows V versus J v , with data points for 
l/kfa = —0.25 (triangle), (circle) and 1.0 (plus sign). Solid 
line shows the GP prediction for a — 1/kf [Eq. (|T0J]. 



tion of 1/kf a. The graph shows that N s increases mono- 
tonically with l/kfa, reflecting that the soliton becomes 
much deeper while its width remains almost the same [f| ■ 
We determine dJ v /dV by measuring V and J v as the 
soliton passes the center of the trap. In the inset in Fig.[3l 
we plot results for l/kfa = —0.25, and 1.0, with a tri- 
angle, circle and plus sign respectively. As expected, the 
result for l/kfa = 1 lies close to the GP prediction [iljj 
(black line) for a = 1/k f and small V, which is 



V 



^TTh 2 n p /4:kfm 2 (n~J v ) . (10) 



Suprisingly, the results for l/kfa = —0.25 and also lie 
near the black line. The variation in dJ v /dV is compara- 
ble to the error in our simulations so, to a good approxi- 
mation, we may take dJ v /dV to be a constant, given by 
Eq. pTJ]) . Given this and mj = 2m, Eq. ([5]) becomes 



(T s /T x ) 



1 



(3/7T 



,1/6 



2/3 



/N s 



ill) 



In Fig. [3] we plot Eq. (fTTj) (solid line) together with 
the numerical results (crosses); the two predictions agree 
within the accuracy of our TDBdG simulations. Figure[3] 
illustrates that the shallow solitons in the BCS regime ac- 
celerate slower than the deep solitons in the BEC regime. 
In the former case, |iV s | is smaller, since more particles 
(unpaired fermions) are present inside the soliton. Ac- 
cording to Eq. (jlip . this implies a larger T s , in agreement 
with our numerical simulations. This is analagous to the 
increase in the period of dark solitons when they are filled 
by an impurity, creating bright-dark solitons 0, Q ■ 

Experimental protocol. We propose an experiment 
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FIG. 4: Proposed experimental scheme to detect the variation 
of T 3 with 1/fc/a. Panels (a) [(c)] and (b) [(d)] show n and 
|A| during protocol A [B]. Horizontal bar shows scale. 

to observe the variation of T s across the BEC-BCS 
crossover. As in Fig.[TJ we consider a. 40 K superfiuid with 
uj x = 27r50 rad s -1 , L± = 3.3 /im and n p = 1.8 x 10 18 
m~ 3 . However, in this case, the initial state does not 
contain a soliton. Instead, we create a soliton dynami- 
cally, using techniques that have been realised in exper- 
iment [9HlH. Firstly, we create a hole in the density of 
the initial state [Fig. Uf a)] by adding a narrow potential 
spike at x — 2.5 /im [20]. Secondly, at t = the cloud at 
x < 2.5 (x > 2.5) /im is imprinted with a phase of zero 
(tt). The potential spike is then smoothly ramped down 
to zero over 0.5 ms (= 0.025T X ) to minimise sound pro- 
duction. We refer to this procedure as protocol A. The 



procedure creates a black soliton plus some sound [20| . 
Since we form the soliton away from the center of the 
trap, it oscillates [Fig. Ufa)]. The amplitude of the order 
parameter is not significantly affected [Fig. HKb)] . 

We now consider a second procedure, referred to as 
protocol B. The soliton is created in the same way but, 
from t = to 10 ms (= 0.5^), we ramp the scatter- 
ing length until 1/fc/fi = —0.5. This reduces the speed 
and depth of the soliton, and Friedel oscillations appear 
in the density profile [Fig. S£c)]. Also, the amplitude 
of the order parameter reduces dramatically [Fig. Hfd)]. 
Unfortunately, the soliton is now too shallow to be ob- 
served experimentally, so we ramp back to 1/fc/a = 1 
from t = 16 (= 0.8T X ) to 26 ms (= 1.3T X ). The soli- 
ton now becomes deeper again, the Friedel oscillations 
disappear [Fig. Hfc)], and the order parameter increases 
[Fig. SJd)]. By comparing the position of the soliton fol- 
lowing protocol B to that following protocol A, the ex- 
perimentalist may prove that T s has been increased. 

Summary. In this work, starting from very general as- 
sumptions, we derive universal relations valid for both 
Bosonic and Fermionic superfluids, relating the soliton 
energy and oscillation period T s to observable quantities. 
Then we solve the TDBdG equations to test these rela- 
tions in a nontrivial theoretical model, and also to pro- 
vide numerical predictions for T s in realistic conditions. 
Our TDBdG simulations, although exceeding challeng- 
ing, represent also a new opportunity to probe unknown 
dynamics across the BEC-BCS crossover. We hope that 
our work, in particular our proposed protocol, will stim- 
ulate experiments to test our predictions. 
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